triangular, backward substitution, crear m, camponormal, las multiplicaciones mtv,mtm



void matriz_esp::_construir_m(matriz<tripla_dbl> & m)
{
        _filas_original = m.cant_filas();
        _columnas_original = m.cant_columnas();
       
        _m_filas = 2 * m.cant_filas() *  m.cant_columnas();
        _m_colum = (m.cant_filas() + 1) * (m.cant_columnas() + 1) - 1;
       
        for(unsigned int i = 0; i < _m_filas; i++)
                _m.push_back(vpid());

        unsigned int fila_i_x = 0;
        unsigned int colu_i_x = 0;
       
        unsigned int fila_i_y = 1;
        unsigned int colu_i_y = 0;

        for(unsigned int i = 0; i < m.cant_filas(); i++)
        {
                for(unsigned int j = 0; j < m.cant_columnas(); j++)
                {
                        //aca los elementos correspondientes a las filas x
                        _m_set_elem(fila_i_x, colu_i_x, -m[i][j].trd());
                        _m_set_elem(fila_i_x, colu_i_x +  m.cant_columnas() + 1, m[i][j].trd());
                                               
                        fila_i_x += 2;
                        colu_i_x++;
       
                        //aca los elementos correspondientes a las filas y
                        _m_set_elem(fila_i_y, colu_i_y, -m[i][j].trd());
                        _m_set_elem(fila_i_y, colu_i_y + 1, m[i][j].trd());
                               
                        fila_i_y += 2;
                        colu_i_y++;
                       
                        //si termine de recorrer una fila en la matriz original debo avanzar una columna mas   
                        if(j ==  m.cant_columnas() - 1)
                        {
                                colu_i_x++;
                                colu_i_y++;
                        }                      
                }
        }
}

void matriz_esp::_construir_mtxm()
{
        /* La matriz resultado tiene la cantidad de filas de Mt por cantidad de columnas de M, que en ambos casos es la cantidad de columnas de M */
       
        _mtxm_filas = _mtxm_columnas = _m_colum;
       
        for(unsigned int i = 0; i < _m_colum; i++)
        {
                _mtxm_pf.push_back(vpid());
                _mtxm_pc.push_back(vpid());
        }
               
        /* Sabemos que la matriz m tiene 2 elementos por cada fila (sino esto no funciona!) */
       
        for(unsigned int i = 0; i < _m_filas; i++)
        {
                /* Estoy en la iesima fila, tengo que hacer 4 multiplicaciones */
               
                vpid fila = _m[i];
                assert(fila.size() == 2);
               
                pair<unsigned int, double> prim = fila[0];
                pair<unsigned int, double> segu = fila[1];
                unsigned int pos_prim = prim.first;
                unsigned int pos_segu = segu.first;
               
                _mtxm_add_elem(pos_prim, pos_prim, prim.second * prim.second);
                _mtxm_add_elem(pos_prim, pos_segu, prim.second * segu.second);
                _mtxm_add_elem(pos_segu, pos_prim, segu.second * prim.second);
                _mtxm_add_elem(pos_segu, pos_segu, segu.second * segu.second);           
        }
}

vector<double> matriz_esp::mult_mt(vector<double> &v)
{
        if(_m_filas != v.size())
                cout << "mult_mt(): Las dimensiones no concuerdan!" << endl;
       
        vector<double> res(_m_colum, 0);
       
        for(unsigned int i = 0; i < _m_filas; i++)
        {
                vpid fila = _m[i];
                assert(fila.size() == 2);
               
                pair<unsigned int, double> prim = fila[0];
                pair<unsigned int, double> segu = fila[1];
                unsigned int pos_prim = prim.first;
                unsigned int pos_segu = segu.first;
               
                res[pos_prim] += v[i] * prim.second;
                res[pos_segu] += v[i] * segu.second;   
        }
       
        return res;    
}


void matriz_esp::_mtxm_triangular()
{
        if(mostrar_debug)
        {
                cout << "Triangulando matriz esparsa..." << endl << endl;
                cout << "  0 %" << endl;
        }
       
        int porcentaje_ant = 0;
       
        /* Version mas copada de gauss by Martin C. y  Martin S. */
               
        for(unsigned int columna_a_anular = 0; columna_a_anular < _m_colum; columna_a_anular++)
        {               //cout << "columna " << columna_a_anular << endl;
                if(mostrar_debug)       // Mensaje de debug, no darle bola
                {
                        double porcentaje = (columna_a_anular * 100) / _m_colum;
                        if(porcentaje > porcentaje_ant)
                        {
                                porcentaje_ant = porcentaje;
                                cout << setw(3);
                                cout << porcentaje << " %" << endl;
                        }
                }              
               
                /* Pongo en 0 la columna a anular y obtengo los multiplicadores */
               
                double elem_diagonal = 0;
               
                for(unsigned int i = 0; i < _mtxm_pc[columna_a_anular].size(); i++)
                {
                        if(_mtxm_pc[columna_a_anular][i].first == columna_a_anular)
                        {
                                elem_diagonal = _mtxm_pc[columna_a_anular][i].second;
                                break;
                        }
                       
                        if(_mtxm_pc[columna_a_anular][i].first > columna_a_anular)
                                break;
                }              
                       
                if(esCero(elem_diagonal))
                {
                        cout << "HAY UN CERO EN LA DIAGONAL " << columna_a_anular << " " << columna_a_anular << endl;
                }
                else
                {      
                        /* El elemento de la diagonal es distinto de 0, creamos los multiplicadores */
                       
                        vector<pid> multiplicadores;

                        for(unsigned int i = _mtxm_pc[columna_a_anular].size() - 1; i > 0; i--)
                        {
                                unsigned int fila = _mtxm_pc[columna_a_anular][i].first;
                               
                                if(fila <= columna_a_anular) break;             /* Si estoy en la diagonal o por arriba, salgo */
                               
                                /* Si sigo por aca es xq estoy por abajo de la diagonal */
                                                               
                                multiplicadores.push_back(pair<unsigned int, double>(fila, (-(_mtxm_pc[columna_a_anular][i].second) / elem_diagonal)));
                               
                                _mtxm_pc[columna_a_anular].pop_back(); // borro el ultimo porque queremos poner ceros en la columna                    
                        }

                        /* Ya tenemos los multipicadores, hay que usarlos. Iteramos las columnas aplicando los multiplicadores */
                       
                        for(unsigned int columna_actual = columna_a_anular + 1; columna_actual < _m_colum ; columna_actual++)
                        {
                                /* Queremos el elemento que esta en la misma fila que el elemento de la diagonal, pero en la columna actual */
                               
                                double constante_multiplicador = 0;
                                                       
                                for(unsigned int i = 0; i < _mtxm_pc[columna_actual].size(); i++)
                                {
                                        if(_mtxm_pc[columna_actual][i].first == columna_a_anular)       //  Encontre el elemento que esta en la columna a modificar y en la fila de la diagonal
                                        {
                                                constante_multiplicador = _mtxm_pc[columna_actual][i].second;
                                                break;
                                        }
                                       
                                        if(_mtxm_pc[columna_actual][i].first > columna_a_anular)        // Me pase y no lo encontre
                                                break;
                                }
                               
                                // Si la constante del multiplicador es 0, significa que las posiciones de la columna actual que estan por debajo de la diagonal no se modifican
                               
                                if(constante_multiplicador == 0)
                                        continue;
                       
                                // Si llegamos hasta aca tenemos que modificar los valores que estan por debajo del elemento de la diagonal en la columna actual
                                                               
                                unsigned int puntero_a_columna = _mtxm_pc[columna_actual].size() - 1;                          
                               
                                for(unsigned int puntero_multiplicadores = 0; puntero_multiplicadores < multiplicadores.size(); puntero_multiplicadores++)
                                {
                                        // El vector de multiplicadores se recorre en orden y la columna al reves
                                       
                                        unsigned int fila_multi = multiplicadores[puntero_multiplicadores].first;
                                        double elemento_actual = 0;                                    
                                       
                                        while(true)
                                        {
                                                if(_mtxm_pc[columna_actual][puntero_a_columna].first == fila_multi)
                                                {
                                                        elemento_actual = _mtxm_pc[columna_actual][puntero_a_columna].second;
                                                        break;                                                 
                                                }
                                               
                                                if(_mtxm_pc[columna_actual][puntero_a_columna].first < fila_multi)
                                                        break;

                                                puntero_a_columna--;
                                        }
                                       
                                        double valor_a_poner = elemento_actual + constante_multiplicador * multiplicadores[puntero_multiplicadores].second;
                                                                               
                                        if(elemento_actual == 0)
                                        {
                                                if(puntero_a_columna == (_mtxm_pc[columna_actual].size() - 1))
                                                {
                                                        _mtxm_pc[columna_actual].push_back(pair<unsigned int,double>(fila_multi,valor_a_poner));
                                                }
                                                else
                                                {
                                                        vector< pid>::iterator it;
                                                        it = _mtxm_pc[columna_actual].begin();
                                                        it = _mtxm_pc[columna_actual].insert(it + puntero_a_columna + 1, pair<unsigned int,double>(fila_multi, valor_a_poner));
                                                }
                                                continue;
                                        }
                                       
                                        if(esCero(valor_a_poner))
                                        {
                                                _mtxm_pc[columna_actual].erase(_mtxm_pc[columna_actual].begin()+ puntero_a_columna );
                                                continue;
                                        }
                                       
                                        /* Si esta todo bien, solo modificamos el valor que existia */
                                                               
                                        _mtxm_pc[columna_actual][puntero_a_columna].second = valor_a_poner;                                    
                                }                                      
                        }      
                       
                        /* Ya modificamos todas las columnas, ahora modificamos el termino independiente */
                       
                        double constante_multiplicador_b = b[columna_a_anular];         // Recordar que columna a anular es la fila donde esta el elemento de la diagonal
                                                                               
                        // Al igual que antes, si la constante del multiplicador es 0, no se cambia nada
                       
                        if(esCero(constante_multiplicador_b))
                                continue;
               
                        // Recorro el vector de multiplicadores, y modifico las posiciones relacionadas con el multiplicador
                       
                        for(unsigned int puntero_multiplicadores = 0; puntero_multiplicadores < multiplicadores.size(); puntero_multiplicadores++)
                        {
                                unsigned int fila_multi = multiplicadores[puntero_multiplicadores].first;                              
                                b[fila_multi] += constante_multiplicador_b * multiplicadores[puntero_multiplicadores].second;
                        }
                }
        }
       
        if(mostrar_debug)
        {
                cout << "100 %" << endl << endl;
                cout << "Matriz esparsa triangulada." << endl;
        }
       
        _mtxm_reestablecer_inv_pf();
}


/* Despeja las incognitas con backward sustitution */

void matriz_esp::_mtxm_backsust()
{
        if(mostrar_debug)
        {
                cout << "Despejando incognitas..." << endl << endl;
                cout << "  0 %" << endl;
        }
       
        int porcentaje_ant = 0;
       
        for(int i = _mtxm_filas - 1; i >= 0; i--)
        {
                if(mostrar_debug)       // Mensaje de debug, no darle bola
                {
                        double porcentaje = ((_mtxm_filas - i - 1) * 100) / _mtxm_filas;
                        if(porcentaje > porcentaje_ant)
                        {
                                porcentaje_ant = porcentaje;
                                cout << setw(3);
                                cout << porcentaje << " %" << endl;
                        }
                }
               
                x[i] = b[i];
               
                int it;
                unsigned int col;              
               
                for(it = _mtxm_pf[i].size() - 1; it >= 0; it--)         // Recorro la fila de atras para adelante
                {                      
                        // Agarro el elemento en la posicion it
                        // Veo en que columna esta el elemento
                        // - Si la columna es mayor que i, procedo
                        // - Si la columna es menor o igual que i, corto
                       
                        if((col = _mtxm_pf[i][it].first) > i)
                                x[i] -= _mtxm_pf[i][it].second * x[col];
                               
                        else
                        {
                                // cout << i << " " << _mtxm_pf[i][it].first << endl;   // Para ver que efectivamente paramos en la diagonal
                                break;
                        }
                }
               
                // Como MtxM no tiene ceros en la diagonal, en _mtxm_pf[i][it] tengo el elemento de la diagonal
                               
                x[i] /= _mtxm_pf[i][it].second;
        }
       
       
       
        if(mostrar_debug)
        {
                cout << "100 %" << endl << endl;
                cout << "Incognitas despejadas." << endl;
        }
}
